#include "headers.hh"
#include "AnaUtils.h"
#include "SysUtils.h"

#define EPSILON 1e-12

Bool_t fkfit = kTRUE;
TH2D *ftpc=0x0;
TH2D *ftof=0x0;
TH2D *fkk=0x0;
TH2D *fzz=0x0;
TH2D *fxi=0x0;
TH1D *fallstat=0x0;
TH2D *fmcmom=0x0;
TH2D *ftpcmc[4]={0x0,0x0,0x0,0x0};
Bool_t fkprint = kFALSE;
Int_t fkHyperMatrix = -999;

const Int_t fcol[]={kBlue, kRed, kGreen+3, kMagenta};
const Int_t fmsty[]={20,21,22,23};

const Int_t fNtype = 4;

Int_t fMatrixYRebin = 1;

TFile * getHyperMatrixFile(const TString pwd, const Int_t khm)
{
  TString tmpdir=pwd;
  TString dirids[100];
  Int_t dircounter=0;
  while(tmpdir.Contains("/")){
    //printf("test0 %s\n", tmpdir.Data());
    tmpdir=tmpdir(tmpdir.First('/')+1,tmpdir.Length());
    //printf("test1 %s\n", tmpdir.Data());
    dirids[dircounter]=tmpdir(0,tmpdir.First('/'));
    tmpdir=tmpdir(tmpdir.First('/'),tmpdir.Length());
    ///printf("test2 %s %d %s\n", tmpdir.Data(), dircounter, dirids[dircounter].Data());
    dircounter++;
  }
  TString d9=dirids[9];
  d9.ReplaceAll("DATA_","");
  d9.ReplaceAll("_coherent","#");
  d9=d9(0,d9.First('#'));
  
  TString stm;
  for(Int_t ii=0; ii<=6; ii++){
    stm+=Form("/%s", dirids[ii].Data());
  }
  stm+=Form("/HyperMatrix/%s/%s.root", d9.Data(), d9.Data());
  
  //===============

  Int_t jpt0=-999;
  Int_t jpt1=-999;

  if(pwd.Contains("JPT5-10")){
    jpt0=5;
    jpt1=10;
  }
  else if(pwd.Contains("JPT10-15")){
    jpt0=10;
    jpt1=15;
  }
  else if(pwd.Contains("JPT15-20")){
    jpt0=15;
    jpt1=20;
  }
  else{
    printf("drawpar pwd bad jpt %s\n", pwd.Data()); exit(1);
  }

  const TString oldjpt=Form("JPT%d-%d", jpt0,  jpt1);
  const TString newjpt=Form("JPT%d-%d", jpt0+1,jpt1);
  stm.ReplaceAll(oldjpt, newjpt);

  TFile *ftm = new TFile(stm);
  return ftm;
}

Double_t Xmin()
{
  return -2;
}

void RebinXall(const Int_t nr)
{
  if(nr==1)
    return;

  TH2D *hhs[]={ftpc, ftof, fkk, fzz, fxi, fmcmom, ftpcmc[0], ftpcmc[1], ftpcmc[2], ftpcmc[3]};
  const Int_t nh=sizeof(hhs)/sizeof(TH2D*);
  for(Int_t ii=0; ii<nh; ii++){
    if(hhs[ii]){
      printf("RebinXall %d: %s %d\n", ii, hhs[ii]->GetName(), nr);
      hhs[ii]->RebinX(nr);
    }
  }

  TH1D *h1s[]={fallstat};
  const Int_t n1=sizeof(h1s)/sizeof(TH2D*);
  for(Int_t ii=0; ii<n1; ii++){
    if(h1s[ii]){
      printf("RebinXall %d: %s %d\n", ii, h1s[ii]->GetName(), nr);
      h1s[ii]->RebinX(nr);
    }
  }  
}

TH2D * PreMatrix(const TString hn, const TString tit, const TH2D * hmat0, const TGraphAsymmErrors * gv1, const Int_t nrby)
{
  TH2D * hmatrix = (TH2D*) hmat0->Clone(hn);
  hmatrix->SetTitle(tit);
  if(nrby>1){
    //hmatrix->RebinY(nrby);

    //rebin and then divide
    for(Int_t ix=1; ix<= hmatrix->GetNbinsX(); ix++){
      for(Int_t iy=1; iy<= hmatrix->GetNbinsY(); iy+=nrby){
        const Double_t cont = hmatrix->Integral(ix,ix, iy, iy+nrby-1)/nrby;
        for(Int_t ib=0; ib<nrby; ib++){
          hmatrix->SetBinContent(ix, iy+ib, cont);
        }
      }
    }

  }
  hmatrix->GetXaxis()->SetRangeUser(gv1->GetX()[0], gv1->GetX()[gv1->GetN()-1]);
  printf("testing ranges %e %e\n", gv1->GetX()[0], gv1->GetX()[gv1->GetN()-1]);
  printf("\n\nprinting hmatrix...\n");
  for(Int_t ix=hmatrix->GetXaxis()->GetFirst(); ix<= hmatrix->GetXaxis()->GetLast(); ix++){
    for(Int_t iy=1; iy<= hmatrix->GetNbinsY(); iy++){
      if(iy>=3 && iy<= hmatrix->GetNbinsY()-1)
        continue;
      printf("%s %s xx %d/%d: %e %e, yy %d/%d: %e\n", hmatrix->GetName(), hmatrix->GetTitle(), ix, hmatrix->GetNbinsX(), hmatrix->GetXaxis()->GetBinLowEdge(ix), hmatrix->GetXaxis()->GetBinUpEdge(ix), iy, hmatrix->GetNbinsY(), hmatrix->GetYaxis()->GetBinCenter(iy)); 
    }
  }
  printf("\n\n");
 
  if(fkprint){
    TCanvas * cc = new TCanvas;
    gPad->SetLogz(0);
    gStyle->SetOptStat(0);

    TH2D * htmp = (TH2D*)hmatrix->Clone("drawingonly");//SysUtils::NormalHist(hmatrix,0,1);

    htmp->Draw("colz");
    cc->Print(hn+Form("_rb%d_khm%d.eps", fMatrixYRebin, fkHyperMatrix));

    SysUtils::ToNaturalScale(htmp);
    htmp->GetXaxis()->SetTitle("p (GeV/c)");
    htmp->Draw("colz");
    cc->Print(hn+Form("_rb%d_khm%d_logx0.eps", fMatrixYRebin, fkHyperMatrix));

    delete cc;
    delete htmp;
  }

  return hmatrix;
}

void CheckTOFstat(const TH1D * all, const TH1D *tpc, TH1D *tof[])
{
  const Int_t np = all->GetNbinsX();
  for(Int_t ii=0; ii<=np+1; ii++){
    Double_t count = 0;
    
    const Double_t nal = all->GetBinContent(ii);
    count += nal;
    
    const Double_t ntp = tpc->GetBinContent(ii);
    count -= ntp;
    
    Double_t nto[]={-999,-999,-999, -999};
    for(Int_t itype=0; itype<fNtype; itype++){
      nto[itype] = tof[itype]->GetBinContent(ii);
      count -= nto[itype];
    }

    printf("CheckTOFstat print ii %d count %.1f -- %s %e %.1f %s %e %.1f %s %e %.1f %.1f %.1f\n", ii, count, all->GetName(), all->GetXaxis()->GetBinCenter(ii), nal, tpc->GetName(), tpc->GetXaxis()->GetBinCenter(ii), ntp, tof[0]->GetName(), tof[0]->GetXaxis()->GetBinCenter(ii), nto[0], nto[1], nto[2]); 
    if(fabs(count)>EPSILON){
      printf("CheckTOFstat Fail!!\n"); exit(1);
    }
  }
}
 
void CorrectTOFfraction(Double_t &p0y, Double_t &e0l, Double_t &e0h,  const Double_t tpcall, const Double_t tofpar, const Double_t tofall)
{
  if(fkprint){
    printf("checkpatch CorrectTOFfraction p0y old %e tpcall %e tofpar %e tofall %e\n", p0y, tpcall, tofpar, tofall);
  }

  p0y = (p0y*tpcall+tofpar)/(tpcall+tofall);

  //have to take into account the correct error propagation, especially when e0lh=0
  e0l = TMath::Sqrt( TMath::Power(e0l*tpcall,2) + tofpar  )/(tpcall+tofall);//correct since tofpar is real count
  e0h = TMath::Sqrt( TMath::Power(e0h*tpcall,2) + tofpar  )/(tpcall+tofall); //correct since tofpar is real count

  if(fkprint){
    printf("checkpatch CorrectTOFfraction p0y new %e tpcall %e tofpar %e tofall %e\n", p0y, tpcall, tofpar, tofall);
  }
}

void CorrectTOFratio(Double_t &p3y, Double_t &e3l, Double_t &e3h, const Double_t tpcpi, const Double_t tofpar,const Double_t  tofpi)
{
  p3y = (p3y*tpcpi+tofpar)/(tpcpi+tofpi);

  //have to take into account the correct error propagation, especially when e0lh=0
  e3l = TMath::Sqrt( TMath::Power(e3l*tpcpi,2) + tofpar )/(tpcpi+tofpi);//correct since tofpar is real count
  e3h = TMath::Sqrt( TMath::Power(e3h*tpcpi,2) + tofpar )/(tpcpi+tofpi);//correct since tofpar is real count
}


//all has the same x-axis
void PatchTOF(TList * ll, TGraphAsymmErrors * tpcFrac[], TGraphAsymmErrors * tpcRatio[], TH1D *tpcCount, TH1D * tofCounts[], TGraphAsymmErrors * tofPIDeff[], const Bool_t kchecksum)
{
  TGraphAsymmErrors * oldsum = new TGraphAsymmErrors(tpcFrac[0]->GetN()); oldsum->SetName(Form("oldsum%s", tpcFrac[0]->GetName())); ll->Add(oldsum);
  TGraphAsymmErrors * newsum = new TGraphAsymmErrors(tpcFrac[0]->GetN()); newsum->SetName(Form("newsum%s", tpcFrac[0]->GetName())); ll->Add(newsum);

  TGraphAsymmErrors * oldfrac[fNtype];
  for(Int_t itype=0; itype<fNtype; itype++){
    oldfrac[itype]=(TGraphAsymmErrors*) tpcFrac[itype]->Clone(Form("tmpoldback%s",tpcFrac[itype]->GetName()));
  }

  for(Int_t itype=0; itype<fNtype; itype++){
    if(fkprint){
      printf("checkpatch PatchTOF itype %d\n", itype);
    }

    if(tofPIDeff){
      tofPIDeff[itype]=new TGraphAsymmErrors;
    }

    if( tpcFrac[itype]->GetN() != tpcRatio[itype]->GetN() ){
      printf("PatchTOF frac and ratio has different N itype %d %d %d\n", itype, tpcFrac[itype]->GetN(), tpcRatio[itype]->GetN()); exit(1);
    }

    for(Int_t ip = 0; ip < tpcFrac[itype]->GetN(); ip++){
      const Double_t xx = tpcFrac[itype]->GetX()[ip];
      if( fabs( xx - tpcRatio[itype]->GetX()[ip] )>EPSILON ){
        printf("PatchTOF frac and ratio has different X itype %d ip %d %e %e\n", itype, ip, tpcFrac[itype]->GetX()[ip], tpcRatio[itype]->GetX()[ip]); exit(1);
      }

      if(fkprint){
        printf("checkpatch ip %d xx %e\n", ip, xx);
      }

      oldsum->SetPoint(ip, xx, oldsum->GetY()[ip]+tpcFrac[itype]->GetY()[ip]);

      const Double_t exl = tpcFrac[itype]->GetEXlow()[ip];
      const Double_t exh = tpcFrac[itype]->GetEXhigh()[ip];

      const Int_t tpcib = tpcCount->GetXaxis()->FindBin(xx);
      if( fabs( xx - tpcCount->GetXaxis()->GetBinCenter(tpcib) ) > EPSILON ){
        printf("PatchTOF tpcCount bin find X wrong itype %d ip %d tpcib %d %e\n",  itype, ip, tpcib, tpcCount->GetXaxis()->GetBinCenter(tpcib)); exit(1); 
      }

      const Double_t tpcall = tpcCount->GetBinContent(tpcib);
      //tpcFrac are updated and can not be used for K/pi calculation because pi is already updated in last itype!!
      const Double_t tpcpi  = tpcall*oldfrac[AnaUtils::kPION]->GetY()[ip]; 

      const Int_t tofib = tpcib; //tofCounts[itype]->GetXaxis()->FindBin(xx);
      if( fabs( xx - tofCounts[itype]->GetXaxis()->GetBinCenter(tofib) ) > EPSILON ){
        printf("PatchTOF tof bin find X wrong itype %d ip %d tofib %d %e\n",  itype, ip, tofib, tofCounts[itype]->GetXaxis()->GetBinCenter(tofib)); exit(1);
      }
      
      const Double_t tofpar = tofCounts[itype]->GetBinContent(tofib);
      const Double_t tofpi  = tofCounts[AnaUtils::kPION]->GetBinContent(tofib);
      Double_t tofall = 0;
      for(Int_t tofii=0; tofii<fNtype; tofii++){
        tofall += tofCounts[tofii]->GetBinContent(tofib);
      }

      //============================> start patching...

      Double_t p0y = tpcFrac[itype]->GetY()[ip]; 
      const Double_t tmp0 = p0y;
      Double_t e0l = tpcFrac[itype]->GetEYlow()[ip];
      Double_t e0h = tpcFrac[itype]->GetEYhigh()[ip];

      Double_t p3y = tpcRatio[itype]->GetY()[ip];
      const Double_t tmp3 = p3y;
      Double_t e3l = tpcRatio[itype]->GetEYlow()[ip];
      Double_t e3h = tpcRatio[itype]->GetEYhigh()[ip];

      if(tofall>EPSILON){        
        CorrectTOFfraction(p0y, e0l, e0h, tpcall, tofpar, tofall);
        tpcFrac[itype]->SetPoint(ip, xx, p0y);
        tpcFrac[itype]->SetPointError(ip, exl, exh, e0l, e0h);
        
        //--
        CorrectTOFratio(   p3y, e3l, e3h, tpcpi,  tofpar, tofpi);
        tpcRatio[itype]->SetPoint(ip, xx, p3y);
        tpcRatio[itype]->SetPointError(ip, exl, exh, e3l, e3h);
      }

      //---
      const Double_t tofeff = (tpcall+tofall)>EPSILON ? tofpar/p0y/(tpcall+tofall) : 0;
      const Double_t toferr = (tpcall+tofall)>EPSILON ? TMath::Sqrt(tofeff*(1-tofeff)/(tofall+tpcall)) : 0;
      if(tofPIDeff){        
        const Int_t tmpn = tofPIDeff[itype]->GetN();
        tofPIDeff[itype]->SetPoint(tmpn, xx, tofeff);
        tofPIDeff[itype]->SetPointError(tmpn, exl, exh, toferr, toferr);
      }
      printf("PatchTOF correcting with tof! itype %d ip %d xx %e mom %e tpcpi %e tpcall %e tofpar %e tofpi %e tofall %e, p0 %e  %e,  p3 %e p3 %e, TOF-PID-EFF %e\n", itype, ip, xx, pow(10,xx),  tpcpi, tpcall, tofpar, tofpi, tofall, tmp0, p0y, tmp3, p3y, tofeff);
      //--
      //<==========================

      newsum->SetPoint(ip, xx, newsum->GetY()[ip]+tpcFrac[itype]->GetY()[ip]);
    }
    TGraphAsymmErrors* oldfr = (TGraphAsymmErrors*) tpcFrac[itype]->Clone(Form("v1%s", tpcFrac[itype]->GetName())); ll->Add(oldfr);   SysUtils::ToNaturalScale(oldfr);
    TGraphAsymmErrors* oldra = (TGraphAsymmErrors*) tpcRatio[itype]->Clone(Form("v1%s", tpcRatio[itype]->GetName())); ll->Add(oldra); SysUtils::ToNaturalScale(oldra);
  }

  SysUtils::ToNaturalScale(oldsum);
  SysUtils::ToNaturalScale(newsum);


  for(Int_t itype=0; itype<fNtype; itype++){
    delete oldfrac[itype];
  }

  //check consistency
  //both loop check numerical identity

  if(kchecksum){
    const Double_t diffeps = EPSILON;

    for(Int_t ii=0; ii<oldsum->GetN(); ii++){
      const Double_t oldii = oldsum->GetY()[ii];
      const Double_t newii = newsum->GetY()[ii];
      
      //counter example:
      //PatchTOF inconsistency in sum of fractions ii 0 oldii 9.949337e-01 newii 9.949357e-01 1.955080e-06
      if(fabs(oldii-newii)>diffeps){
        printf("PatchTOF inconsistency in sum of fractions ii %d oldii %e newii %e %e\n", ii, oldii, newii, newii-oldii); exit(1);
      }
    }
    
    for(Int_t itype=0; itype<fNtype; itype++){
      for(Int_t ii=0; ii<tpcRatio[itype]->GetN(); ii++){
        const Double_t fi  = tpcFrac[itype]->GetY()[ii];
        const Double_t fpi = tpcFrac[AnaUtils::kPION]->GetY()[ii];
        if(fpi<EPSILON){
          continue;
        }
        const Double_t rr = fi/fpi;
        const Double_t cr = tpcRatio[itype]->GetY()[ii];
        if(fabs(rr-cr)>diffeps){
          printf("PatchTOF inconsistency in to-pion ratio! itype %d ii %d fi %e fpi %e rr %e cr %e %e\n", itype, ii, fi, fpi, rr, cr, rr-cr); exit(1);
        }
      }
    }
  }
}

void PrintTOF(/*TGraphAsymmErrors * frac[], */ TGraphAsymmErrors *eff[], TGraphAsymmErrors *gp0[], const Bool_t kmctr)
{
  TCanvas * cc = new TCanvas;
  gStyle->SetOptStat(0);
  gPad->SetLogx(1);

  /*
  for(Int_t itype=0; itype<fNtype; itype++){
    frac[itype]->SetMinimum(0);
    frac[itype]->SetMaximum(0.6);
    frac[itype]->GetXaxis()->SetLimits(0.1,10);
    frac[itype]->GetXaxis()->SetTitle("p (GeV/c)");
    frac[itype]->GetYaxis()->SetTitle("TOF fraction");
    frac[itype]->SetTitle("");
    frac[itype]->SetLineColor(fcol[itype]);
    frac[itype]->SetMarkerStyle(fmsty[itype]);
    frac[itype]->SetMarkerColor(fcol[itype]);

    frac[itype]->Draw(itype?"lp":"alp");
  }
  gPad->SetLogy(1);
  cc->Print(Form("TOFfracMCtr%d.eps", kmctr));
  */

  gPad->SetLogy(0);
  for(Int_t itype=0; itype<fNtype; itype++){
    eff[itype]->SetMinimum(0);
    eff[itype]->SetMaximum(0.7);
    eff[itype]->GetXaxis()->SetLimits(0.1,10);
    eff[itype]->GetXaxis()->SetTitle("p (GeV/c)");
    eff[itype]->GetYaxis()->SetTitle("TOF PID efficiency");
    eff[itype]->SetTitle("");
    eff[itype]->SetLineColor(fcol[itype]);
    eff[itype]->SetMarkerStyle(fmsty[itype]);
    eff[itype]->SetMarkerColor(fcol[itype]);

    eff[itype]->Draw(itype?"lp":"alp");
  }
  cc->Print(Form("TOFeffMCtr%d_khm%d.eps", kmctr, fkHyperMatrix));

  TList *tmplist = new TList;
  tmplist->SetOwner(kTRUE);

  for(Int_t itype=0; itype<fNtype; itype++){
    TGraphAsymmErrors * upg= (TGraphAsymmErrors*) gp0[itype]->Clone(Form("upg%d", itype)); tmplist->Add(upg);
    SysUtils::ToNaturalScale(upg);
    //TGraphAsymmErrors * old= (TGraphAsymmErrors*) oldp0[itype]->Clone(Form("old%d", itype)); tmplist->Add(old);
    //SysUtils::ToNaturalScale(old);
    upg->SetMinimum(0);
    upg->SetMaximum(1);
    upg->GetXaxis()->SetLimits(0.1, 10);
    upg->GetXaxis()->SetTitle("p (GeV/c)");
    upg->GetYaxis()->SetTitle("uncorrected fraction");
    upg->SetLineColor(fcol[itype]);    
    upg->SetMarkerStyle(fmsty[itype]);
    upg->SetMarkerColor(fcol[itype]);

    upg->SetLineStyle(kSolid);
    upg->SetTitle("black dotted: #slash{TOF} vs. color solid: +TOF");

    //old->SetMarkerStyle(fmsty[itype]);
    //old->SetLineStyle(kDotted);

    upg->Draw(itype?"lp":"alp");
    //old->Draw("lp");
  }
  cc->Print(Form("comFracMCtr%d_khm%d.eps",kmctr, fkHyperMatrix));

  delete tmplist;

  gPad->SetLogx(0);
  delete cc;
}

TGraphAsymmErrors * getDelta(const TGraphAsymmErrors * grdata, const TGraphAsymmErrors * grmctruth)
{
  const TGraphAsymmErrors *grs[]={grdata, grmctruth};
  Int_t starts[]={-999,-999};
  Int_t stops[]={-999,-999};
  const Int_t npt = SysUtils::GetCommonLimits(2, grs, starts, stops);

  TGraphAsymmErrors * grdelta = new TGraphAsymmErrors;
  for(Int_t ii=0; ii<npt; ii++){
    const Int_t idata = starts[0]+ii;
    const Int_t imct  = starts[1]+ii;

    const Double_t xdata = grdata->GetX()[idata];
    if(fabs(xdata-grmctruth->GetX()[imct])>EPSILON){
      printf("drawpar TGraphAsymmErrors * getDelta xdata and xmctruth not match!! ii %d idata %d imct %d -- xdata %e xmct %e\n", ii, idata, imct, xdata, grmctruth->GetX()[imct]); exit(1);
    }

    const Double_t ydelta = grdata->GetY()[idata] - grmctruth->GetY()[imct];
    grdelta->SetPoint(ii, xdata, ydelta);

    const Double_t xlerr = grdata->GetEXlow()[idata];
    const Double_t xherr = grdata->GetEXhigh()[idata];
    const Double_t ylerr = grdata->GetEYlow()[idata];
    const Double_t yherr = grdata->GetEYhigh()[idata];
    grdelta->SetPointError(ii, xlerr, xherr, ylerr, yherr);
  }

  return grdelta;
}

void IniTOF(TH1D * &htpc, TH1D * htofcount[], TList *ll)
{
  htpc = ftpc->ProjectionX("tpcstat",0, ftpc->GetNbinsY()+1); 
      
  TH1D *htoffrac[]={0x0,0x0,0x0,0x0};
      
  TGraphAsymmErrors * gtoffrac[]={0x0,0x0,0x0,0x0};

  for(Int_t itype=0; itype<fNtype; itype++){
    const Int_t tmpy = ftof->GetYaxis()->FindBin(itype+EPSILON);
    htofcount[itype]=ftof->ProjectionX(Form("tofCount%d", itype), tmpy, tmpy); 
    
    htoffrac[itype]=(TH1D*) htofcount[itype]->Clone(Form("tofFrac%d", itype)); ll->Add(htoffrac[itype]);
    htoffrac[itype]->Divide(htoffrac[itype], fallstat, 1, 1, "B");
        
    gtoffrac[itype]=SysUtils::HistToGraph(Form("toffrac%d", itype), htoffrac[itype]); ll->Add(gtoffrac[itype]); SysUtils::ToNaturalScale(gtoffrac[itype]);
  }
      
  CheckTOFstat(fallstat, htpc, htofcount);
}

void FitGraph(TList * ll, const TString tn, TGraphAsymmErrors *gp00[], TGraphAsymmErrors *gp01[], TGraphAsymmErrors *gp02[], TGraphAsymmErrors *gp03[])
{
  //chi2
  if(tn.Contains("fith")){
    TGraphAsymmErrors *grchi2=SysUtils::TreeToGraph(tn,"1",(Char_t*)"xx",(Char_t*)"chidof"); ll->Add(grchi2);
  }
  
  //totn from horiginalyield
  TGraphAsymmErrors *gtotn=SysUtils::TreeToGraph(tn,"1",(Char_t*)"xx",(Char_t*)"totn"); ll->Add(gtotn);
  
  //all log10p base
  for(Int_t itype=0; itype<fNtype; itype++){
    //mn2d->Scan("le00:he00:ea0:le00-he00:le00-ea0") 
    //old: all equal
    //new: le00=he00, !=ea0 -> ea0 is set from Binomial stat.

    //gp00[itype]=SysUtils::TreeToGraph(tn,Form("xx>%e", Xmin()),(Char_t*)"xx",Form("p%d0", itype),   (Char_t*)"",(Char_t*)"",Form("le%d0", itype),Form("he%d0", itype),0); gp00[itype]->SetName(Form("gr%sp%d00", tn.Data(), itype)); ll->Add(gp00[itype]);
    gp00[itype]=SysUtils::TreeToGraph(tn,Form("xx>%e", Xmin()),(Char_t*)"xx",Form("p%d0", itype),   (Char_t*)"",(Char_t*)"",Form("ea%d", itype),Form("ea%d", itype),0); gp00[itype]->SetName(Form("gr%sp%d00", tn.Data(), itype)); ll->Add(gp00[itype]);

    //test stat error  -->
    for(Int_t itmpgr=0; itmpgr<gp00[itype]->GetN(); itmpgr++){
      const Double_t tot =  gtotn->GetY()[itmpgr];
      const Double_t frac = gp00[itype]->GetY()[itmpgr];
      const Double_t stat =gp00[itype]->GetEYlow()[itmpgr];
      const Double_t expected = sqrt(frac*(1-frac)/(tot+EPSILON));
      printf("teststat tot %e frac %e stat %e exp %e delta %e\n", tot, frac, stat, expected, stat-expected);
      //gp00[itype]->SetPointError(itmpgr, gp00[itype]->GetEXlow()[itmpgr], gp00[itype]->GetEXhigh()[itmpgr], expected, expected);
    }
    //<---

    gp01[itype]=SysUtils::TreeToGraph(tn,Form("xx>%e", Xmin()),(Char_t*)"xx",Form("p%d1", itype),   (Char_t*)"",(Char_t*)"",Form("le%d1", itype),Form("he%d1", itype),0); gp01[itype]->SetName(Form("gr%sp%d01", tn.Data(), itype)); ll->Add(gp01[itype]);
    SysUtils::ToNaturalScale(gp01[itype]);
    gp02[itype]=SysUtils::TreeToGraph(tn,Form("xx>%e", Xmin()),(Char_t*)"xx",Form("p%d2", itype),   (Char_t*)"",(Char_t*)"",Form("le%d2", itype),Form("he%d2", itype),0); gp02[itype]->SetName(Form("gr%sp%d02", tn.Data(), itype)); ll->Add(gp02[itype]);
    SysUtils::ToNaturalScale(gp02[itype]);
    gp03[itype]=SysUtils::TreeToGraph(tn,Form("xx>%e", Xmin()),(Char_t*)"xx",Form("rr%d", itype),   (Char_t*)"",(Char_t*)"",Form("er%d", itype), Form("er%d", itype),0);  gp03[itype]->SetName(Form("gr%sp%d03", tn.Data(), itype)); ll->Add(gp03[itype]);
 
    gp00[itype]->SetLineColor(fcol[itype]); gp00[itype]->SetMarkerStyle(fmsty[itype]); gp00[itype]->SetMarkerColor(fcol[itype]);
    gp01[itype]->SetLineColor(fcol[itype]); gp01[itype]->SetMarkerStyle(fmsty[itype]); gp01[itype]->SetMarkerColor(fcol[itype]);
    gp02[itype]->SetLineColor(fcol[itype]); gp02[itype]->SetMarkerStyle(fmsty[itype]); gp02[itype]->SetMarkerColor(fcol[itype]);
    gp03[itype]->SetLineColor(fcol[itype]); gp03[itype]->SetMarkerStyle(fmsty[itype]); gp03[itype]->SetMarkerColor(fcol[itype]);
  }

  //print gp00: p0: xx=log10P
  printf("\n\nprinting p10, i.e. p0:p, just for type0...\n");
  for(Int_t ii=0; ii<gp00[0]->GetN(); ii++){
    printf("gp00[0] ii %d/%d log10P %e %e %e, p0 %e %e %e\n", ii, gp00[0]->GetN(), gp00[0]->GetX()[ii], gp00[0]->GetEXlow()[ii], gp00[0]->GetEXhigh()[ii], gp00[0]->GetY()[ii], gp00[0]->GetEYlow()[ii], gp00[0]->GetEYhigh()[ii]); 
  }
  printf("\n\n");

  //before TOF correction, used in PrintTOF, no NaturalScale yet
  for(Int_t itype=0; itype<fNtype; itype++){
    TGraphAsymmErrors* old00 = (TGraphAsymmErrors*) gp00[itype]->Clone(Form("v0%s", gp00[itype]->GetName())); ll->Add(old00); SysUtils::ToNaturalScale(old00);
    TGraphAsymmErrors* old03 = (TGraphAsymmErrors*) gp03[itype]->Clone(Form("v0%s", gp03[itype]->GetName())); ll->Add(old03); SysUtils::ToNaturalScale(old03);
  }
}

void MCTruthGraph(TList*ll, TGraphAsymmErrors * gp50[], TGraphAsymmErrors * gp51[], TGraphAsymmErrors * gp52[], TGraphAsymmErrors * gp53[], const TGraphAsymmErrors* gfit)
{
  Int_t parbin[fNtype];
  for(Int_t itype=0; itype<fNtype; itype++){
    parbin[itype] = fmcmom->GetYaxis()->FindBin(itype+EPSILON);
    printf("MCTruthGraph parbin %d: %d\n", itype, parbin[itype]);
  }
  TH1D * hmcall = fmcmom->ProjectionX(Form("%sall", fmcmom->GetName()), parbin[0], parbin[3]);
  if(hmcall->GetSumw2N()==0){
    hmcall->Sumw2();
  }
  TH1D * hmcpi  = fmcmom->ProjectionX(Form("%spi",  fmcmom->GetName()), parbin[AnaUtils::kPION], parbin[AnaUtils::kPION]);
  if(hmcpi->GetSumw2N()==0){
    hmcpi->Sumw2();
  }
  ll->Add(hmcall);
  ll->Add(hmcpi);

  Double_t xmin = -999;
  Double_t xmax = -999;
  if(gfit){
    xmin = gfit->GetX()[0]-EPSILON;
    xmax = gfit->GetX()[gfit->GetN()-1]+EPSILON;
    printf("MCTruthGraph xmin xmax %e %e\n", xmin, xmax);
  }
  else{
    for(Int_t ib=1; ib<=hmcpi->GetNbinsX(); ib++){
      if(hmcpi->GetBinContent(ib)>EPSILON && xmin<-990){
        xmin = hmcpi->GetBinCenter(ib);
        break;
      }
      }
    for(Int_t ib=hmcpi->GetNbinsX(); ib>=1; ib--){
      if(hmcpi->GetBinContent(ib)>EPSILON && xmax<-990){
        xmax = hmcpi->GetBinCenter(ib);
        break;
      }
    }
  }
  printf("MCTruthGraph xmin %e xmax %e\n", xmin, xmax);

  for(Int_t itype=0; itype<fNtype; itype++){
    TH1D * htype50 = fmcmom->ProjectionX(Form("%sp%d50", fmcmom->GetName(), itype), parbin[itype], parbin[itype]);
    ll->Add(htype50);
    if(htype50->GetSumw2N()==0){
      htype50->Sumw2();
    }
    htype50->Divide(htype50, hmcall, 1, 1, "B");
    gp50[itype] = SysUtils::HistToGraph(Form("grmn2dp%d%d", itype, 50), htype50, -1, 0x0, xmin, xmax); 
    ll->Add(gp50[itype]);
    //save before any modification
    TGraphAsymmErrors* old50 = (TGraphAsymmErrors *) gp50[itype]->Clone(Form("v0%s", gp50[itype]->GetName())); ll->Add(old50); SysUtils::ToNaturalScale(old50);

    TH1D * htype53 = fmcmom->ProjectionX(Form("%sp%d53", fmcmom->GetName(), itype), parbin[itype], parbin[itype]);
    ll->Add(htype53);
    //htype53->Divide(htype53, hmcpi,  1, 1, "B");
    //same as getMCETAall.C
    for(Int_t ibin=1; ibin<=htype53->GetNbinsX(); ibin++){
      const Double_t ydown = hmcpi->GetBinContent(ibin);
      const Double_t yup = htype53->GetBinContent(ibin);

      if(ydown<EPSILON || yup < EPSILON){
        htype53->SetBinContent(ibin,0); htype53->SetBinError(ibin,0);
        continue;
      }
      const Double_t rat = yup/ydown;
      const Double_t err = rat * TMath::Sqrt(1/yup+1/ydown); //correct since yup and ydown are real counts
      htype53->SetBinContent(ibin, rat); htype53->SetBinError(ibin,err);
      printf("testtest %i %e %e %e %e\n", ibin, yup, ydown, rat, err);
    }
    //<--

    gp53[itype] = SysUtils::HistToGraph(Form("grmn2dp%d%d", itype, 53), htype53, -1, 0x0, xmin, xmax); 
    ll->Add(gp53[itype]);
    TGraphAsymmErrors* old53 = (TGraphAsymmErrors *) gp53[itype]->Clone(Form("v0%s", gp53[itype]->GetName())); ll->Add(old53); SysUtils::ToNaturalScale(old53);

    if(gp50[itype]->GetN()!=gp53[itype]->GetN()){
      printf("MCTruthGraph gp50[itype]->GetN()!=gp53[itype]->GetN() %d %d\n",gp50[itype]->GetN(), gp53[itype]->GetN()); 
      for(Int_t itt=0; itt<gp50[itype]->GetN(); itt++){
        printf("50 itype %d %d: %e %e\n", itype, itt, gp50[itype]->GetX()[itt], gp50[itype]->GetY()[itt]);
      }
      printf("\n");
      for(Int_t itt=0; itt<gp53[itype]->GetN(); itt++){
        printf("53 itype %d %d: %e %e\n", itype, itt, gp53[itype]->GetX()[itt], gp53[itype]->GetY()[itt]);
      }
      exit(1);
    }

    gp50[itype]->SetLineColor(fcol[itype]); gp50[itype]->SetMarkerStyle(fmsty[itype]); gp50[itype]->SetMarkerColor(fcol[itype]);
    gp53[itype]->SetLineColor(fcol[itype]); gp53[itype]->SetMarkerStyle(fmsty[itype]); gp53[itype]->SetMarkerColor(fcol[itype]);

    if(ftpcmc[itype]){
      TH1D * hamp = 0x0; TH1D * hmpv = 0x0; TH1D * hwid = 0x0; TH1D * hchi = 0x0; TH1D * hres = 0x0;
      SysUtils::FitSlicesY(ftpcmc[itype],  hamp, hmpv, hwid, hres, hchi, 1, 0,1,0x0); 
      //gp50[itype] = SysUtils::HistToGraph(Form("fitgrmn2dp%d%d", itype, 50), hamp, -1, 0x0, xmin, xmax); ll->Add(gp50); SysUtils::ToNaturalScale(gp50);
      gp51[itype] = SysUtils::HistToGraph(Form("grmn2dp%d%d", itype, 51), hmpv, 20, 0x0, xmin, xmax); ll->Add(gp51[itype]); SysUtils::ToNaturalScale(gp51[itype]);
      gp52[itype] = SysUtils::HistToGraph(Form("grmn2dp%d%d", itype, 52), hres, -1, 0x0, xmin, xmax); ll->Add(gp52[itype]); SysUtils::ToNaturalScale(gp52[itype]); 
      delete hamp;
      delete hmpv;
      delete hres;
      delete hwid;
      delete hchi; 
      gp51[itype]->SetLineColor(fcol[itype]); gp51[itype]->SetMarkerStyle(fmsty[itype]); gp51[itype]->SetMarkerColor(fcol[itype]);
      gp52[itype]->SetLineColor(fcol[itype]); gp52[itype]->SetMarkerStyle(fmsty[itype]); gp52[itype]->SetMarkerColor(fcol[itype]);      
    }

  }
}

//log10p base!
void EtaCorrection(TList *ll, TGraphAsymmErrors * gr0s[], TGraphAsymmErrors * gr3s[])
{
  //to-pion ratio after correction is totally off!!
  return;

  /*
  //tuned for IDFF, maxeta = 1e-6 is the best between 0.9, 0.6, 0.4, 0.2, 0.01. As maxeta decreases, the consistency between etas converges to the direct truth
  //1e-6; has bias, eta mean is lower than eta integrated
  const TString pwd = gSystem->pwd();
  Double_t maxeta = -999;
  //both use 0.9, the nominal value for track eta cut
  if(pwd.Contains("THnIncl")){
    maxeta = 0.9;
  }
  else if(pwd.Contains("THnIDFF")){
    maxeta = 0.9;
  }
  else{
    printf("EtaCorrection wrong pwd %s\n", pwd.Data()); exit(1);
  }
  printf("EtaCorrection maxeta %e\n", maxeta);

  for(Int_t itype=0; itype<fNtype; itype++){
    if(gr0s[itype]->GetN()!=gr3s[itype]->GetN()){
      printf("EtaCorrection itype %d n0!=n3 %d %d\n", itype, gr0s[itype]->GetN(), gr3s[itype]->GetN()); exit(1);
    }

    for(Int_t ip=0; ip<gr0s[itype]->GetN(); ip++){
      const Double_t xx = gr0s[itype]->GetX()[ip];
      if(TMath::Abs(xx-gr3s[itype]->GetX()[ip])>EPSILON){
        printf("EtaCorrection gr0 and gr3 have different x! itype %d ip %d xx %e %e\n", itype, ip, xx, gr3s[itype]->GetX()[ip]); exit(1);
      }
      const Double_t mom = TMath::Power(10, xx);
      const Double_t eta = 0.5*(AnaUtils::Eta0()+AnaUtils::Eta1());
      const Double_t pt = SysUtils::ToPt(mom, eta);

      const Double_t fcorr = AnaUtils::FractionEtaCorrection(itype, pt, maxeta);
      SysUtils::ScaleGraphPoint(gr0s[itype], ip, fcorr);

      const Double_t rcorr = AnaUtils::RatioEtaCorrection(itype, pt);
      SysUtils::ScaleGraphPoint(gr3s[itype], ip, rcorr);

      //printf("EtaCorrection itype %d ip %d xx %e mom %e eta %e pt %e fcorr %e rcorr %e\n", itype, ip, xx, mom, eta, pt, fcorr, rcorr);
    }
    TGraphAsymmErrors * old0 =(TGraphAsymmErrors*) gr0s[itype]->Clone(Form("v2%s",gr0s[itype]->GetName())); ll->Add(old0); SysUtils::ToNaturalScale(old0);
    TGraphAsymmErrors * old3 =(TGraphAsymmErrors*) gr3s[itype]->Clone(Form("v2%s",gr3s[itype]->GetName())); ll->Add(old3); SysUtils::ToNaturalScale(old3);
  }
  */
}

void recalToPionRatio(TGraphAsymmErrors * frac, TGraphAsymmErrors *pionfr, TGraphAsymmErrors *topi, TList *ll)
{
  //modify y value only, keep the error by weighted mean

  TGraphAsymmErrors * oldtopi = (TGraphAsymmErrors *) topi->Clone(Form("%swm", topi->GetName()));
  ll->Add(oldtopi);

  for(Int_t ip=0; ip<topi->GetN(); ip++){
    const Double_t xx = topi->GetX()[ip];

    if( fabs(frac->GetX()[ip]-xx)>EPSILON || fabs(pionfr->GetX()[ip]-xx)>EPSILON ){
      printf("recalToPionRatio xx error ip %d xx %e %e %e\n", ip, xx, frac->GetX()[ip], pionfr->GetX()[ip]); exit(1);
    }

    const Double_t newyy = frac->GetY()[ip]/pionfr->GetY()[ip];

    topi->SetPoint(ip, xx, newyy);

    oldtopi->SetPoint(ip, xx, newyy/oldtopi->GetY()[ip]);
  }
}

void drawpar(const TString tn, const TString outname)
{
  TList *ll = new TList;

  ll->Add(ftpc);
  if(fkk) ll->Add(fkk);
  if(fzz) ll->Add(fzz);
  if(fxi) ll->Add(fxi);
  if(ftof) ll->Add(ftof);
  if(fmcmom) ll->Add(fmcmom);


  //====================================================================================================
  //set up TOF  
  TH1D * htpc = 0x0;
  TH1D * htofcount[]={0x0,0x0,0x0,0x0};

  if(ftof){
    //check TOFstat, get TOF fraction!
    IniTOF(htpc, htofcount, ll);
  }

  //========================================================================
  //gpXY: Form("gr%sp%dXY", tn.Data(), itype)
  //
  //X: x-axis
  //[0-4 data]
  //0: p
  //1: pT (weighted mean)
  //2: z (weighted mean)
  //3: xi (weighted mean)
  //
  //[5-9 MC true]
  //5: MC true p base
  //6: MC true pT base
  //7: MC true z base
  //8: MC true xi base
  //
  //4: p base data - MC truth, or data/truth
  //9: pT base data - MC truth, or data/truth
  //
  //Y: y-axis
  //0: uncorrected fraction
  //1: mean
  //2: sigma
  //3: to-pion ratio
  //4: cov wrt pion
  //8: TOF eff
  //========================================================================

  //fit results
  //p0, p1, p2, p3 =rr(ratio), p4 =cova(cov wrt pion) vs xx
  TGraphAsymmErrors *gp00[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors *gp01[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors *gp02[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors *gp03[]={0x0,0x0,0x0,0x0};

  Int_t filebad = 0;
  
  if(fkfit){
    FitGraph(ll, tn, gp00, gp01, gp02, gp03);
    
    //============================================
    //correct with TOF from sources gp00, gp03!! all has x-axis log10P, so NO natural scale to hists
    if(ftof){
      //------------------
      //correct gp00, gp03
      TGraphAsymmErrors * gtofeff[]={0x0,0x0,0x0,0x0};
      PatchTOF(ll, gp00, gp03, htpc, htofcount, gtofeff, kFALSE);
      
      for(Int_t itype=0; itype<fNtype; itype++){
        ll->Add(gtofeff[itype]);
        gtofeff[itype]->SetTitle(gtofeff[itype]->GetName());
        gtofeff[itype]->SetName(Form("gr%sp%d08", tn.Data(), itype));
        SysUtils::ToNaturalScale(gtofeff[itype]);
      }
      
      if(fkprint){
        PrintTOF(gtofeff, gp00, kFALSE);
      }
    }

    // //check file ok after tof patching  --------------->
    // for(Int_t ipoint=0; ipoint < gp00[0]->GetN(); ipoint++){
    //   //pion must be largest
    //   if(gp00[1]->GetY()[ipoint]<gp00[0]->GetY()[ipoint] || gp00[1]->GetY()[ipoint]<gp00[2]->GetY()[ipoint] || gp00[1]->GetY()[ipoint]<gp00[3]->GetY()[ipoint]){
    //     printf("drawpar pwd %s file not OK failtype1 ipoint %d frac %e %e %e %e\n", gSystem->pwd(), ipoint, gp00[0]->GetY()[ipoint], gp00[1]->GetY()[ipoint], gp00[2]->GetY()[ipoint], gp00[3]->GetY()[ipoint]);                                          
    //     filebad = 1;
    //   }
    //   //depend on the separation of p/K fractions which is different in data and MC
    //   /*
    //   const Double_t pthres = TMath::Log10(1);
    //   if(gp00[2]->GetX()[ipoint] > pthres && gp00[2]->GetY()[ipoint] < gp00[0]->GetY()[ipoint]){
    //     printf("drawpar pwd %s file not OK failtype2 ipoint %d frac %e %e %e %e\n", gSystem->pwd(), ipoint, gp00[0]->GetY()[ipoint], gp00[1]->GetY()[ipoint], gp00[2]->GetY()[ipoint], gp00[3]->GetY()[ipoint]);
    //     filebad = 2;
    //   }
    //   */
    // }
    
    // //kaon > proton for last point, given that last point is > 15 GeV
    // const Double_t pthres = TMath::Log10(15);
    // const Int_t ilast = gp00[2]->GetN()-1;
    // if(gp00[2]->GetX()[ilast] > pthres){
    //   if(gp00[2]->GetY()[ilast] < gp00[0]->GetY()[ilast]){
    //     printf("drawpar pwd %s file not OK failtype2 ilast %d frac %e %e %e %e\n", gSystem->pwd(), ilast, gp00[0]->GetY()[ilast], gp00[1]->GetY()[ilast], gp00[2]->GetY()[ilast], gp00[3]->GetY()[ilast]);
    //     filebad = 2;
    //   }
    // }
    
    // //<---------------------

    //EtaCorrection(ll, gp00, gp03);
  }

  //=========================================================================================================================================

  //MC true 5: mom base
  TGraphAsymmErrors * gp50[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * gp51[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * gp52[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * gp53[]={0x0,0x0,0x0,0x0};

  if(fmcmom){
    MCTruthGraph(ll, gp50, gp51, gp52, gp53, gp00[0]);

    if(ftof){
      TGraphAsymmErrors * mctofeff[]={0x0,0x0,0x0,0x0};
      //correct 50, 53
      PatchTOF(ll, gp50, gp53, htpc, htofcount, mctofeff, kTRUE);

      for(Int_t itype=0; itype<fNtype; itype++){
        ll->Add(mctofeff[itype]);
        mctofeff[itype]->SetTitle(mctofeff[itype]->GetName());
        mctofeff[itype]->SetName(Form("gr%sp%d58", tn.Data(), itype));
        SysUtils::ToNaturalScale(mctofeff[itype]);
      }

      if(fkprint){
        PrintTOF(mctofeff, gp50, kTRUE);
      }
    }

    //print 50
    for(Int_t tmpit=0; tmpit<fNtype; tmpit++){
      for(Int_t ip=0; ip<gp50[tmpit]->GetN(); ip++){
        printf("drawpar type50: tmpit %d ip %d/%d %e %e\n", tmpit, ip, gp50[tmpit]->GetN(), gp50[tmpit]->GetX()[ip], gp50[tmpit]->GetY()[ip]);
      }
      //print 53
      for(Int_t ip=0; ip<gp53[tmpit]->GetN(); ip++){
        printf("drawpar type53: tmpit %d ip %d/%d %e %e\n", tmpit, ip, gp53[tmpit]->GetN(), gp53[tmpit]->GetX()[ip], gp53[tmpit]->GetY()[ip]);
      }
    }
    //EtaCorrection(ll, gp50, gp53);
  }

  //=========================================================================================================================================

  //print matrixPt x: log10P, y: log10Pt
  if(fkk && fzz && fxi){

    //pt can not be rebinned because it is needed for outlier removal. p based can not be used there because of different x bins!
    TH2D * fmatrixPt = PreMatrix("matrixLog10Pt", "; log10 p (GeV/c) ; log10 p_{T} (GeV/c)",           fkk, fkfit?gp00[0]:gp50[0], 1); ll->Add(fmatrixPt);
    
    //print matrixZ x: log10P, y: z
    TH2D * fmatrixZ  = PreMatrix("matrixZ",       "; log10 p (GeV/c) ; z = track_{pT}/jet_{pT}",       fzz, fkfit?gp00[0]:gp50[0], fMatrixYRebin); ll->Add(fmatrixZ);
    
    //print matrixXi x: log10P, y: xi
    TH2D * fmatrixXi = PreMatrix("matrixXi",      "; log10 p (GeV/c) ; #xi = ln(jet_{pT}/track_{pT})", fxi, fkfit?gp00[0]:gp50[0], fMatrixYRebin); ll->Add(fmatrixXi);
    
    Int_t npoint = -999;
    const Double_t thres = 50;

    const Bool_t kPtAcceptance = kFALSE;    
    const Bool_t kZXiAcceptance = kFALSE;

    //{1,2,3,6,7,8}3
    TGraphAsymmErrors * ktopi[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
    TGraphAsymmErrors * ptopi[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
    TGraphAsymmErrors * etopi[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};

    //{1,2,3,6,7,8}0
    TGraphAsymmErrors * pifrac[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
    TGraphAsymmErrors * kfrac[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
    TGraphAsymmErrors * pfrac[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
    TGraphAsymmErrors * efrac[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0};

    for(Int_t itype=0; itype<fNtype; itype++){
      //no cut edge for p->pT
      TGraphAsymmErrors * gp10=0x0;
      TGraphAsymmErrors * gp13=0x0;
      if(fkfit){
        gp10 = SysUtils::MeasurementConversion(Form("gr%sp%d10", tn.Data(), itype), gp00[itype], fmatrixPt, kPtAcceptance, thres); if(!gp10){printf("gp10 null!!\n"); exit(1);} 
        ll->Add(gp10);
        SysUtils::ToNaturalScale(gp10);
        if(npoint<0){
          npoint = gp10->GetN();
        }
        else if(npoint!=gp10->GetN()){
          printf("drawpar npoint not consistent itype %d npoint %d GetN %d\n", itype, npoint, gp10->GetN()); exit(1);
        }
        
        gp13 = SysUtils::MeasurementConversion(Form("gr%sp%d13", tn.Data(), itype), gp03[itype], fmatrixPt, kPtAcceptance, thres); if(!gp13){printf("gp13 null!!\n"); exit(1);}
        ll->Add(gp13);
        SysUtils::ToNaturalScale(gp13);
        
        TGraphAsymmErrors * gp20 = SysUtils::MeasurementConversion(Form("gr%sp%d20", tn.Data(), itype), gp00[itype], fmatrixZ, kZXiAcceptance, thres); 
        ll->Add(gp20);
        
        TGraphAsymmErrors * gp23 = SysUtils::MeasurementConversion(Form("gr%sp%d23", tn.Data(), itype), gp03[itype], fmatrixZ, kZXiAcceptance, thres); 
        ll->Add(gp23);
        
        TGraphAsymmErrors * gp30 = SysUtils::MeasurementConversion(Form("gr%sp%d30", tn.Data(), itype), gp00[itype], fmatrixXi, kZXiAcceptance, thres); 
        ll->Add(gp30);
        
        TGraphAsymmErrors * gp33 = SysUtils::MeasurementConversion(Form("gr%sp%d33", tn.Data(), itype), gp03[itype], fmatrixXi, kZXiAcceptance, thres); 
        ll->Add(gp33);

        if(itype==1){
          pifrac[0]=gp10;
          pifrac[1]=gp20;
          pifrac[2]=gp30;
        }
        else if(itype==2){
          kfrac[0]=gp10;
          kfrac[1]=gp20;
          kfrac[2]=gp30;
        }
        else if(itype==0){
          pfrac[0]=gp10;
          pfrac[1]=gp20;
          pfrac[2]=gp30;
        }
        else if(itype==3){
          efrac[0]=gp10;
          efrac[1]=gp20;
          efrac[2]=gp30;
        }

        if(itype==2){
          ktopi[0]=gp13;
          ktopi[1]=gp23;
          ktopi[2]=gp33;
        }
        else if(itype==0){
          ptopi[0]=gp13;
          ptopi[1]=gp23;
          ptopi[2]=gp33;
        }
        else if(itype==3){
          etopi[0]=gp13;
          etopi[1]=gp23;
          etopi[2]=gp33;
        }
      }
      
      //MC
      if(fmcmom){
        TGraphAsymmErrors * gp60 = SysUtils::MeasurementConversion(Form("gr%sp%d60", tn.Data(), itype), gp50[itype], fmatrixPt, kPtAcceptance, thres); if(!gp60){printf("gp60 null!!\n"); exit(1);}
        ll->Add(gp60);
        SysUtils::ToNaturalScale(gp60);
        
        TGraphAsymmErrors * gp63 = SysUtils::MeasurementConversion(Form("gr%sp%d63", tn.Data(), itype), gp53[itype], fmatrixPt, kPtAcceptance, thres); if(!gp63){printf("gp63 null!!\n"); exit(1);}
        ll->Add(gp63);
        SysUtils::ToNaturalScale(gp63);
        
        TGraphAsymmErrors * gp70 = SysUtils::MeasurementConversion(Form("gr%sp%d70", tn.Data(), itype), gp50[itype], fmatrixZ, kZXiAcceptance, thres);
        ll->Add(gp70);
        
        TGraphAsymmErrors * gp73 = SysUtils::MeasurementConversion(Form("gr%sp%d73", tn.Data(), itype), gp53[itype], fmatrixZ, kZXiAcceptance, thres);
        ll->Add(gp73);
        
        TGraphAsymmErrors * gp80 = SysUtils::MeasurementConversion(Form("gr%sp%d80", tn.Data(), itype), gp50[itype], fmatrixXi, kZXiAcceptance, thres);
        ll->Add(gp80);
        
        TGraphAsymmErrors * gp83 = SysUtils::MeasurementConversion(Form("gr%sp%d83", tn.Data(), itype), gp53[itype], fmatrixXi, kZXiAcceptance, thres);
        ll->Add(gp83);
        
        if(itype==1){
          pifrac[3]=gp60;
          pifrac[4]=gp70;
          pifrac[5]=gp80;
        }
        else if(itype==2){
          kfrac[3]=gp60;
          kfrac[4]=gp70;
          kfrac[5]=gp80;
        }
        else if(itype==0){
          pfrac[3]=gp60;
          pfrac[4]=gp70;
          pfrac[5]=gp80;
        }
        else if(itype==3){
          efrac[3]=gp60;
          efrac[4]=gp70;
          efrac[5]=gp80;
        }

        if(itype==2){
          ktopi[3]=gp63;
          ktopi[4]=gp73;
          ktopi[5]=gp83;
        }
        else if(itype==0){
          ptopi[3]=gp63;
          ptopi[4]=gp73;
          ptopi[5]=gp83;
        }
        else if(itype==3){
          etopi[3]=gp63;
          etopi[4]=gp73;
          etopi[5]=gp83;
        }

        //don't use, depending on X3 which is not finial here
        /*
        if(fkfit){
          //-->
          //get 40 = 00-50, 43 = 03-53
          TGraphAsymmErrors * gp40 = getDelta(gp00[itype], gp50[itype]); gp40->SetName(Form("gr%sp%d40", tn.Data(), itype)); ll->Add(gp40);
          TGraphAsymmErrors * gp43 = getDelta(gp03[itype], gp53[itype]); gp43->SetName(Form("gr%sp%d43", tn.Data(), itype)); ll->Add(gp43); 
          
          //get 90 = 10-60, 93 = 13-63
          TGraphAsymmErrors * gp90 = getDelta(gp10, gp60); gp90->SetName(Form("gr%sp%d90", tn.Data(), itype)); ll->Add(gp90);
          TGraphAsymmErrors * gp93 = getDelta(gp13, gp63); gp93->SetName(Form("gr%sp%d93", tn.Data(), itype)); ll->Add(gp93);
          
          TGraphAsymmErrors * gp41=0x0;
          TGraphAsymmErrors * gp42=0x0;
          TGraphAsymmErrors * dummy=0x0;//nothing should happen do dummy
          if(SysUtils::SimpleScaleToRef(gp01[itype], 0x0, gp51[itype], gp41, dummy, kFALSE)){
            gp41->SetName(Form("gr%sp%d41", tn.Data(), itype)); ll->Add(gp41); gp41->SetLineColor(fcol[itype]); gp41->SetMarkerStyle(fmsty[itype]); gp41->SetMarkerColor(fcol[itype]);
          }
          if(SysUtils::SimpleScaleToRef(gp02[itype], 0x0, gp52[itype], gp42, dummy, kFALSE)){
            gp42->SetName(Form("gr%sp%d42", tn.Data(), itype)); ll->Add(gp42); gp42->SetLineColor(fcol[itype]); gp42->SetMarkerStyle(fmsty[itype]); gp42->SetMarkerColor(fcol[itype]);
          }
          //<---
        }
        */
      }
    }

    const Int_t nrecal = 6;
    for(Int_t irecal=0; irecal<nrecal; irecal++){
      //void recalToPionRatio(TGraphAsymmErrors * frac, TGraphAsymmErrors *pionfr, TGraphAsymmErrors *topi, TList *ll)
      if(!kfrac[irecal])
        continue;

      recalToPionRatio(kfrac[irecal], pifrac[irecal], ktopi[irecal], ll);
      recalToPionRatio(pfrac[irecal], pifrac[irecal], ptopi[irecal], ll);
      recalToPionRatio(efrac[irecal], pifrac[irecal], etopi[irecal], ll);
    }
  }

  //after PrintTOF because gp00 and oldgp00 have both to be in log10
  for(Int_t itype=0; itype<fNtype; itype++){
    if(fkfit){
      SysUtils::ToNaturalScale(gp00[itype]);
      SysUtils::ToNaturalScale(gp03[itype]);
    }

    if(fmcmom){
      SysUtils::ToNaturalScale(gp50[itype]);
      SysUtils::ToNaturalScale(gp53[itype]);
    }
  }
  
  TFile *fout=new TFile(Form("%spar%s%s_khm%d.root", filebad? Form("badfile%d", filebad):"", outname.Data(), Form("_rb%d", fMatrixYRebin), fkHyperMatrix),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc == 5){
    fkprint = atoi(argv[4]);
  }
  else if(argc != 4){
    printf("drawpar argc != 4\n"); 
    return 1;
  }

  fMatrixYRebin=atoi(argv[2]);
  fkHyperMatrix = atoi(argv[3]);

  const TString pwd = gSystem->pwd();

  //=====

  //original.root
  TFile *fhistroot = new TFile("../../hist.root");
  if(!fhistroot->IsOpen()){
    printf("no ../hist.root!!!!\n");
    exit(1);
  }
 
  TH2D *zzdumy = (TH2D*)fhistroot->Get("zz");
  if(!zzdumy){
    printf("drawpar error in original.root no zzdumy!\n"); exit(1);
  }
  fallstat =  zzdumy->ProjectionX("allstat", 0, zzdumy->GetNbinsY()+1);  

  ftpc = (TH2D*)fhistroot->Get("horiginal");
  if(!ftpc){
    printf("drawpar error in original.root no horiginal!\n"); fhistroot->ls(); exit(1);
  }

  ftof=0x0;
  if(!pwd.Contains("TOF0")){
    ftof = (TH2D*)fhistroot->Get("tof");
    if(ftof){
      printf("drawpar using tof!!\n");
    }
  }

  fmcmom =(TH2D*) fhistroot->Get("mcmom");
  if(fmcmom){
    printf("drawpar fmcmom!\n");
  }

  for(Int_t itype=0; itype<fNtype; itype++){
    ftpcmc[itype]=(TH2D*) fhistroot->Get(Form("tpcmc%d", itype));
    if(!ftpcmc[itype]){
      printf("no tmpcmc %d\n", itype); //exit(1);
    }
  }

  //=====
  TFile * fHyperMatrix = 0x0;
  if(fkHyperMatrix==0){
    fHyperMatrix = fhistroot;
  }
  else{
    fHyperMatrix  = getHyperMatrixFile(pwd, fkHyperMatrix);
    if(!fHyperMatrix->IsOpen()){
      printf("fHyperMatrix not open!! %d %s %s\n", fkHyperMatrix, pwd.Data(), fHyperMatrix->GetName()); exit(1);
    }
  }
  printf("pleasecheck HyperMatrix %s\n", fHyperMatrix->GetName());

  fkk = (TH2D*)fHyperMatrix->Get("kk");// pt:tpcip
  if(!fkk){
    printf("drawpar error in original.root no kk!\n"); //fHyperMatrix->ls(); exit(1);
  }

  fzz = (TH2D*)fHyperMatrix->Get("zz");// pt:tpcip
  if(!fzz){
    printf("drawpar error in original.root no zz!\n"); //fHyperMatrix->ls(); exit(1);
  }

  fxi = (TH2D*)fHyperMatrix->Get("xi");// pt:tpcip
  if(!fxi){
    printf("drawpar error in original.root no xi!\n"); //fHyperMatrix->ls(); exit(1);
  }


  Int_t xrebin = 1;
  if(pwd.Contains("XB9"))
    xrebin = 9;
  else if(pwd.Contains("XB7"))
    xrebin = 7;
  else if(pwd.Contains("XB5"))
    xrebin = 5;
  else if(pwd.Contains("XB3"))
    xrebin = 3;
  else{
    printf("drawpar wrong XB! %s\n", pwd.Data()); exit(1);
  }
  printf("drawpar xrebin %d\n", xrebin);

  RebinXall(xrebin);

  //===================================================

  const TString fin(argv[1]);
  TFile *ffit = TFile::Open(fin);
  fkfit = kTRUE;
  if(!ffit || !ffit->IsOpen()){
    fkfit = kFALSE;
  }
  printf("drawpar main fkfit %d\n", fkfit);

  const TString tn = fin.Contains("mn2d")?"mn2d":"fith";

  TString outname(fin);
  outname.ReplaceAll("/","");
  outname.ReplaceAll(".","");
  outname.ReplaceAll("out","");
  outname.ReplaceAll("_post","");
  outname.ReplaceAll("_save","");
  outname.ReplaceAll(".root","");
  outname.ReplaceAll("*","");

  drawpar(tn, outname);

  printf("drawpardone!!\n");
  return 0;
}
